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Abstract. We investigate the propagation of wave-packets on graphene in a 
perpendicular magnetic field and the appearance of collapses and revivals in the time- 
evolution of an initially localised wave-packet. The wave-packet evolution in graphene 
differs drastically from the one in an electron gas and shows a rich revival structure 
similar to the dynamics of highly excited Rydberg states. We present a novel numerical 
wave-packet propagation scheme in order to solve the effective single-particle Dirac- 
Hamiltonian of graphene and show how the collapse and revival dynamics is affected 
by the presence of disorder. Our effective numerical method is of general interest for 
the solution of the Dirac equation in the presence of potentials and magnetic fields. 
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1. Introduction 

Since its experimental discovery, graphene has become a hot topic in sohd state 
physics [H [2]. On the theoretical side, the similarities between the effective single- 
particle Hamiltonian for graphene in the low-energy regime and the Dirac equation for 
massless particles have been noted [31 HJ [5] and lead to the prediction of "relativistic 
effects" in graphene. One example is the non-transient zitterbewegung in the presence 
of an external magnetic field [HI E]- In order to study the zitterbewegung and other 
transport effects in graphene one needs to solve the time-dependent Dirac equation 
with high accuracy and efficiency. Previous time-dependent propagation methods have 
often been restricted to eigenstate decompositions. However, the evaluation of integrals 
involving a finite set of eigenf unctions is cumbersome and also limited to energetically 
low lying eigenstates. Additionally the eigenstate approach can only be used if the 
analytic solutions are known, since the numerical determination of a full set of energy 
eigenstates is in general not feasible. 

In section [2] we show that the prior knowledge of the eigenstates is not required in 
order to obtain a reliable solution of the time-dependent Dirac-Hamiltonian of graphene. 
Our algorithm is accurate up to the machine error and stable for any propagation time. 
In addition it is very flexible and thus allows us to calculate the time-evolution in 
inhomogeneous magnetic fields and arbitrary potential landscapes. In section [3] we 
study quantitatively the scattering in realistic potentials and obtain the fidelity (which 
is related to the mobility) of graphene in a magnetic field. For wave-packets with large 
initial momenta we reveal a rich collapse and revival structure, which has a strong 
similarity to the dynamics of highly excited states of Rydberg atoms [H [9] . We derive 
analytic results for the semiclassical expectation value of the position operator and the 
autocorrelation function which are in very good agreement with the numerical results. 
Subsequently, we study the effect of an impurity potential and analyse the suppression 
of the autocorrelation function corresponding to a rapid fidelity decay. Interestingly, 
revivals in the centre-of-mass coordinate survive the disorder-perturbations if the 
correlation length of the potential is on the order of the cyclotron radii of the occupied 
states. In general, our algorithm provides a basic method to model time-dependent 
transport through mesoscopic systems [IDl [HI [12] • 

2. Wave-packet propagation for the graphene-system 

The electronic structure of graphene is obtained from the tight binding model of the 
two-dimensional honeycomb lattice. In the resulting band structure two inequivalent 
vanishing bandgaps (often referred to as different valleys or Dirac points) appear, which 
are located at the corners of the Brillouin zone [131 [T3]. The band structure around each 
Dirac point is conveniently described in the framework of an effective k ■ p perturbation 
theory and the resulting Hamiltonian [15] 
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has the same structure as the 2 + 1 dimensional Dirac Hamiltonian, where the velocity 
is given by v ~ lO^m/s, and a^, ay denote the corresponding Pauli matrices |16j . 
The linear approximation is valid for energies close to the Dirac points. For higher 
energies the trigonal warping gets important [17]. We neglect this effect here, since we 
only consider wave packets which can be decomposed into eigenstates from — 0.7eV to 
0.7 eV. The Hamiltonians ([1]) at the two Dirac points are obtained by switching the 
sign of r = ±1. Here, we study the bulk region of a sheet of graphene and therefore we 
can neglect effects of the boundaries which may couple the valleys and induce symmetry 
breaking [18l[l9]. 

In the following we present a new algorithm for the propagation of a wave-packet on 
a sheet of graphene in an anisotropic magnetic field B-^{r) oriented perpendicularly to 
the sheet. The Hamiltonian can also contain a position-dependent potential landscape 
V{r). The resulting system is described by the following Hamiltonian 

^ / \^(r) —ivhkx + vhky — veA{r) 

\ —ivhkx — vhky — veA*{r) y(r) 

For simplicity we only use the Hamiltonian with r = 1, but it is straightforward to apply 
the same algorithm to the second Dirac point. The vector potential of the magnetic 
field A{r) = {Ax{r), Ay{r),0) is connected to the magnetic field B{r) = V x A{r) and 
enters the Hamiltonian as A{r) = Ax{r) — iAy{r). As initial state of the system we 
choose an arbitrary two spinor wavefunction ip{t = 0) = {ipA{t = 0),4'Bit = 0)}. The 
action of the time-evolution operator on the wavefunction is given by 

|^(t))=exp(-^m)|V;(0)). (3) 

In principle, there exist several possibilities to calculate the action of the time-evolution 
operator. For a system with a discrete set of eigenstates defined by 7i|0„) = Sn\4>n) a 
favourable way is given by the expansion in eigenstates. In this case, the time-evolution 
becomes 

= T.i'Pnim) exp (-^^nt) (4) 

Such a propagation algorithm was already used in references [6l [7] to evaluate the 
zitterbewegung in graphene and also for other semiconductor materials [20]. Due to the 
lack of analytical known solutions of the Dirac equation, this scheme can only be applied 
to a limited set of problems. Thus it is favourable to develop an algorithm which does 
not require any knowledge of the eigenstates and can be applied to arbitrary mesoscopic 
systems. A similar starting point exists in quantum chemistry, where no analytic 
solutions are known for the molecular potential energy surfaces. In all these cases it 
is favourable to use an algorithm which directly gives the time-evolution of the system 
without evaluating all eigenstates first [21]. All relevant information about stationary 
properties, i.e. the energy spectrum or the Green's function and the local density of 
states, is encoded in time-dependent observables, which are extracted afterwards from 
the time-evolved wave-packet [22l [TO]. 
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For the time-dependent solution of the Dirac Hamiltonian, we use a rehable and 
efficient algorithm based upon the polynomial expansion of the propagator in Chebyshev 
polynomials [23]. The main idea of the algorithm is to expand the time-evolution 
operator on the interval [—1,1] using a set of Chebyshev polynomials, which are 
particularly suitable since they distribute the numerical error of the expansion equally 
on the whole interval. As a first step the eigenvalues of the Hamiltonian 7i have to be 
projected into the interval of the expansion by introducing a normalised Hamiltonian 

^ T~C /I 

The scaling-energy AE has to be big enough to cover all eigenenergies contributing to 
the initial wave-packet. With the aid of ([5]) the time-evolution is given by 

m)) = E «n T„(-ii7no,^) 1^(0)) , (6) 

n=0 \ '^il' / " V ' 

where T„ stands for the series of Chebyshev polynomials of the ffist kind. The expansion 
coefficients are given by 

with the Bessel function J„. Only the action of the polynomial of Hnovm on the initial 
state is needed to calculate the time-evolution. Thus it is sufficient to recursively 
calculate the vectors P„ by the following modified recursion relations for Chebyshev 
polynomials: 

n = 1^(0)), (8) 

Pi = 4.orm|^(0)), (9) 
Pn= — 2iifnorm-Pn-l + Pn~2- (10) 

In order to avoid the fermion doubling problem of a simple ffist-order nearest-neighbour 
derivative [23] we use a Fourier representation of the spinors and apply the kinetic energy 
operator in momentum space [21]. Another possibility is to introduce certain non-local 
ffist-order derivatives [25]. In our case, the action of the normalised Hamiltonian acting 
on a two spinor wavefunction (V^a, V'b) is given by 

^(^)^A - V (hj^-^ {ih -ky)J' + ei(r)) 
~ Ae[ V{r)iPB - V {hj^~' {ik, + ky)J^ + ei*(r)) V^a 

where and JF stand for the two-dimensional Fast Fourier Transform. This type of 
propagation algorithm allows us to evaluate the time-evolution of the wave-packet for 
arbitrary long propagation times, as long as the wavefunction does not leave the region 
of interest. A good measure of the propagation quality is given by the comparison 
of the analytic autocorrelation function for the homogeneous magnetic field with the 
numerical result. Our results show that the maximum error is indeed on the order 
of the unavoidable finite machine precision. In contrast to all other propagation 
algorithms for graphene we are able to incorporate arbitrarily shaped potentials and 
also inhomogeneous magnetic fields in our configurations. 
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3. Wave-packet revivals 
3.1. Model and observables 

In the following we apply the propagation algorithm to an infinite sheet of graphene 
put in a perpendicular magnetic field B{r) = {0,0, B) and set the potential V{r) to 
zero. This system was proposed by Rusin et. al. [6j and Schliemann [7] to study the 
non-transient zitterbewegung of a wave-packet on graphene. The eigenstates of the Dirac 
Hamiltonian in a perpendicular magnetic field are known and given by the harmonic 
oscillator eigenfunctions [26]. In the lowest Landau level, the ground state wavefunction 
has Gaussian shape and is located on one sublattice. As initial condition we consider a 
kicked Gaussian of the form 

The width is taken as ao = h/{eB) and the initial momentum k can point along any 
direction in the x — y plane, but is set to fc = koBy in the following. The pseudospin ak 
of the wave-packet is determined by the phase-relation between the two spinor entries. 
For (p = 71 /2, the wave-packet is mostly electron-like since the momentum is pointing 
along the y-direction. Consequently the wave-packet is mostly hole-like ii ip = 37r/2 and 
a mixture of both states for intermediate angles. 

During the propagation we track several observables, including the expectation 
value of the position operator 

r{t) = m)\r\m)- (12) 
We refer to r{t) as the "centre-of-mass" although the calculated motion describes 
massless particles. The trembling motion of the centre-of-mass is visible in the videos 
linked from figure [T](a) for different initial setups. Another important property, the 
local pseudospin, is given by {ip{r,t)\(Tk\ip{r,t)) at each gridpoint r and colour coded 
in the video. The local pseudospin reveals the local electron or hole character of the 
wave-packet and shows an anti-symmetry with respect to the x-axis if the pseudospin 
of the initial state is oriented orthogonal to the initial momentum. An advantage of 
our method is the extraction of stationary properties, like the local density of states 
(LDOS) and the Green's function, out of the time-evolution [10]. For example the 
relative strength of each energy eigenstate with respect to the initial state is given by 
the Fourier transform of the autocorrelation function 

C{t) = {^|J{t = 0)\m)■ (13) 
The same concept allows us to calculate the E x B drift motion and eigenstates in the 
quantum Hall regime of graphene for high bias currents [26] . The LDOS is given by the 
Laplace transformation of the autocorrelation function 

n{E) oc -53 / dt e'^*/^C(t). (14) 



As shown in figure [Tj^b) the LDOS of the trembling wave-packet reveals pronounced 
peaks at the energies of the Landau levels £n = v sign(n) w2e-B^|n|. The peak height 
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Figure 1. (a) Some frames of the time-evolution of a Icicked Gaussian wave-packet 
for different initial pscudospin configurations {B = 5T, ko = a/ eB/Ti, oq = h/{eB)). 
The red pin stands for the centre-of-mass, the height for the absolute value and the 
colour for the local pseudospin given by (■)/'('", t)\ak\ijj{r, t)). A video of the propagation 
is available from http: //www. quantumdynamics.de/graphene for initial pseudospin 
ak 11 k (moviela.mp4, 3.0MB) and ak ± k (movielb.mp4, 3.2MB). (b) Local density 
of states calculated from the wave-packet propagation. The spectrum is antisymmetric 
if the initial pseudospin is parallel to the momentum and symmetric if pscudospin and 
momentum are orthogonal to each other. Red dots denote the analytical energies of 
the Landau levels and coincide up to a very high accuracy with the numerical positions. 



(which is proportional to the enclosed area in the LDOS) indicates the overlap of the 
initial Gaussian wave-packet with the different Hermite polynomials of the Landau 
levels. For (p = 7r/2, mostly electrons are occupied whereas for ip = n holes and electrons 
are equally occupied. Our computational scheme derives the overlap information 
without calculating the numerically unstable integrals over Hermite poljTiomials used 
in previous methods. Therefore we can increase the initial momentum fco of the wave- 
packet and study the quantum mechanical propagation of high Landau levels. 

The trembling motion on graphene is present since both, electron and hole-like 
states, contribute to the initial wave-packet. For high initial momenta the index of the 
average Landau level is given by 

= 4*- <i5) 

which follows from the semiclassical quantisation condition flA.16|) . Our quantum 
mechanical calculations support this assertion as shown in the LDOS of figure [21 A 
strongly kicked Gaussian wave-packet mostly occupies one energy branch if the initial 
pseudospin is parallel to the initial momentum. For an initial momentum which puts 
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Figure 2 

Landau level of n 



Spectrum of a kicked Gaussian wave-packet with an average occupied 
24.5 (B 



IT, ko = 7y/eB/h). The hole-like contribut ions are 
approximately 400 times smaller than the electron-like parts. The red dots show the 
position of the analytical Landau levels as well as the probability distribution p7|) used 
for the following approximations. 



the Gaussian wave-packet into the 24.5* Landau level, the hole contributions are 400 
times smaller compared to the electron parts. In such a scenario the trembling motion 
may be neglected and other interesting revival phenomena appear. In atomic physics 
revival effects are commonly observed for valence electrons in the Coulomb potential of 
Rydberg atoms and have been studied in theory [8] as well as in experiments [H [27] . 
Reference [28] contains a recent review about revivals of quantum wave-packets. In order 
to get an idea about the rich revival structure of a quantum wave-packet on graphene, 
it is instructing to look at the numerical results in figure [3] and the corresponding 
video. Wave-packet revivals are typically not present in semiconducting material in a 
perpendicular magnetic field since there the dynamics is governed by a purely quadratic 
Hamiltonian. Graphene and toplogical insulators form their own class due to their 
non-quadratic Hamilton operators. 



3.2. Revivals in the centre- of -mass position 

As first phenomena we study the collapse and the revival of the centre-of-mass of a 
graphene wave-packet. The quantum mechanical position of the centre-of-mass cannot 
be found analytically. Thus we introduce a classical picture where the centre-of-mass 
is given by a weighted sum over the centre-of-masses moving along quantised cyclotron 
orbits. The centre-of-mass of the classical cyclotron motion in the n-th Landau level is 
defined by z{t) = x(t) + iy{t) with 

Zn{t) = Le'''-' (16) 

as derived in Appendix A The occupation P„ depends on the initial momentum ko, the 



initial pseudospin and the width ao of the wave-packet. As an approximation we find 
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Figure 3. Revivals and fractional revivals of a cyclotron wave-packet in a magnetic 
field of -B = 1 T. The red pin stands for the quantum mechanical centre-of-mass, 
whereas the green pin is calculated with the semiclassical expression p^ . The 
initial wave-packet in (a) was chosen such that the contributing eigenstates follow 
the Gaussian distribution pT]) with n = 60, cr = 1, leading to a Poincare cyclotron 
time of to = Tc\. During the time-evolution the centre-of-mass collapses to the centre of 
a circular wavefunction and later shows pronounced revivals. For example (b) displays 
a quarter revival at t = l/2nTci, (c) a third revival at t = 2/3nTc\ and (d) a half 
revival at t = nTd- All revivals have specific symmetries which also occur in the 
subsequent revivals. After a certain time a localised wave-packet emerges in (e) as a 
mirror revival at t = 2nTci. In this case, the quantum mechanical centre-of-mass is at 
the opposite side compared to the semiclassical centre-of-mass. This phenomenon has 
already been seen in Rydberg revivals [29] . (f ) The full revival where the semiclassical 
and the quantum mechanical centre-of-mass coincide occurs for t = AnTd- The 
related movie (movie2 .mp4, 3.0 MB) which illustrates the dynamics is available from 
http : / / www . quEintumdynamics . de /graphene 



from numerical calculations for high average Landau levels a Gaussian distribution of 
the form 

= exp ^ , (17) 



shown in figure [2l Here n stands for the average occupied Landau level and a denotes the 
Landau level spread. The sum over the Landau levels is normalised if n ^ a is fulfilled. 
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Figure 4. Collapse of a Dirac wave-packet (n = 60, tr = 1, B = IT). Both panels 
show the time-evolution of the x- and the y-position of the wave packet. The black 
dots are obtained from a quantum mechanical calculation and fit quite well to the 
classical centre-of-mass which is shown by the red line and given by the semiclassical 
approximation ()20|1 . A small deviation is visible which stems from the fact that the 
centre-of-mass of the initial wave-packet has already a strong angular dispersion (see 
figure [3|) and thus the average cyclotron radius / is overestimated. 



In this limit the cyclotron radii /„ do not differ substantially from the cyclotron radius 
of the average occupied Landau level / := In. Thus, the time-evolution of the average 
position is approximated by 



iu)„t 



Z{t)=l J2 



The last sum contains both, the centre-of-mass collapse as well as the revivals. Numerical 
results of this classical centre-of-mass motion are shown in figure [3l In the corresponding 
video the classical approximation is marked by the green pin, which almost perfectly 
coincides with the exact quantum mechanical result indicated by the red pin. 



3.2.1. Collapse: Next, we extract the times of the initial collapse of the centre-of-mass 
and the revivals out of equation (ITSl) . We expand the cyclotron frequency uOn around 
the average Landau level fi: 



u)[n 




(19) 



For high Landau levels the energy spectrum becomes denser and the sum can be 
approximated by an integral which we evaluate analytically after the series is reduced 
to the first two terms: 

z{t) = l / rfnP(n)e'^(")* 
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I , ( (n-nV\ / eB 1 / In-n 
dnexp — — — exp w\ ——r= 1 — - 



= / exp iiunt) ■ exp I r I • (20) 

Equation (120!) shows that the centre-of-mass motion completes one cyclotron orbit 
during the cyclotron time of the average Landau level fi 



271 2h, , 
V \l en 

and also a Gaussian decay on the timescale 



4 / h 

Tcoii = — (22) 
va V eB 

Our numerical calculations in figure H] show the same behaviour. Very small deviations 
from the exact result stem from the initial radial dispersion of the Gaussian wave-packet, 
which is not included in the classical theory. 

3.2.2. Revivals: The revival of the classical centre-of-mass takes place if all particles 
orbiting along the contributing Landau orbits reunite at the initial position. We 
calculate the revival time from the Taylor expansion of the classical centre-of-mass 
motion inserted in equation ( fTSi) . In order to meet the conditions for a revival 
\z{T^cv(cu))\ = /, all elements of the sum 

oo 

^(^) ~ ^ 51 ^-^P (i"^"^ + K'^ ~ ^)^'n^ + ...) (23) 

n=— oo 

must have the same phase. Thus, the n-dependent addend of the exponent inuj'^t must 
be a multiple of 27ri. This leads to a revival time of the classical centre-of-mass given by 



^ 27r 47r /2/i,_,, , , 

-^rev(CM) = —r = — \ = 2 n Tel. (24) 

oj'^ V ^ eB 

In this case the classical centre-of-mass of the wave-packet is on the opposite 
side compared to the centre-of-mass of the quantum mechanical calculation. This 
phenomenon is visible in the comparison between the classically and the quantum 
mechanically calculated expectation value of the x coordinate for the first revival time 
shown in figure [5]^a). Exactly the same effect occurs in the mirror revivals of Rydberg 
atoms [291 [30]. the second revival the classical and the quantum mechanical centre- 
of-mass coincide again as is evident from figure [5](b). 

3. 3. Revivals in the autocorrelation function 

In contrast to the centre-of-mass revivals, the revivals in the autocorrelation function can 
be understood without a semiclassical approximation. By definition the time-evolution 
of the autocorrelation function is given by 

oo 

C{t) = {(Pit = O)|0(t)) = e-'^"*/'*P(n). (25) 
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Figure 5. (a) First and (b) second revival of the centre-of-mass (n = 60, cr = 1, 
B = IT). Red line shows the centre-of-mass given by the semiclassical approximation 
()18|) . This result is compared with a full quantum mechanical calculation shown by 
the black dots. Both results fit well with the remarkable difference for odd revivals. 
The quantum mechanical result is opposing the classical result, which has already been 
seen in Rydberg atoms. 



Again, we expand the square root dependence in a Taylor series 

£{n) = vpeBh\n\ 1 - ' + ' 

* V 2 n 8 16 n'^ 



leading to an approximate autocorrelation function 



oo 

C{t)= ^-i{enHn-n)S'^^{n-nfS'i+l(n-n)^S'i' + ...)t/hp^ 



n] 



(26) 



(27) 



where stands for the first derivative of the dispersion and £'^, denote the higher 
order derivatives. It is astonishing that all orders are prominently encoded in the 
quantum mechanical motion and influence the picture at the relevant revival times. 
From the term i£nt/h in the exponent of equation (1271) we obtain the phase oscillation 
of the autocorrelation function on the timescale 
2'iTn 



T 

-J- rtf 



2h 



28 



\£n\ V y eB\n\ 

This oscillation takes place when the wave-packet moves over the initial position. At 
the classical cyclotron time 

27rn 27r [2h\n\ 



\£' 



eB 



(29) 



the wave-packet comes back to the initial position and thus the autocorrelation function 
approximately reaches its initial value, as shown in figure [6]^a). For a revival the terms 
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Figure 6. Revivals in the autocorrelation function of a graphene cyclotron wave- 
packet (n = 20 a = 1, B = IT), (a) The real part shows fast oscillations on the 
timcscalc Tqsc which exhibit revivals every Td corresponding to the time of a classical 
cyclotron orbit. During the time in panel (b) the dispersion of the initially localised 
wave-packet grows, leading to a uniform distribution which recovers to a full revival 
after T^q^ • 



proportional to the second derivative have to be a multiple of 2tt. This leads to a revival 
for 



27rn Stt 2h\n\^ 
T.e. = 2^ = -^-LL = 4|n|T., (30) 

which is best seen in the time-evolution of the absolute value of the autocorrelation 
function (figure EJ^b)). Of course, the revival hierarchy continues to higher orders. In 
atomic physics the next time is the so called "super revival time" 



27rh IGtt 2h\n\^ _o , , 

-/super = 6t-— 7 = \ — = 8n Tcl = 2\n\ Trcv (31) 

\c^'\ V \ eB 

The huge range of timescales requires a highly accurate wave-packet propagation, since 
the phase of the propagated wave-packet needs to be accurate for at least 4|np phase 
oscillations in the autocorrelation function. The used polynomial propagation is well 
suited for this task since the accumulated error is reduced by using only very few but 
long timesteps. 

3.4- Fractional revivals 

Interestingly, we can also distinguish fractional revivals occurring between the full 
revivals [31] which are again encoded in the autocorrelation function [32]. Fractional 
revivals have been experimentally measured in Rydberg atoms for the first time [33]. The 
underlying phase effect can be also employed for prime number factorisation [3ll [35] . 
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Figure 7. Poincare sections of tlie autocorrelation function in order to study fractional 
revivals {n = 60, a = 1, B = T). The autocorrelation function is calculated for 
successive classical cyclotron times n and suitable subdivisions p/q as shown in (a), 
(b) The resulting autocorrelation function shows full revivals, a mirror revival and 
several fractional revivals. There the absolute value features local extremal points 
which are labelled in the plot. 



Recently several experiments have implemented this scheme into a NMR setup |36j . 
cold atoms [37] or optically [38] and successfully factorised numbers. Note that the 
phase properties for wave-packets in graphene differ from those of the aforementioned 
experiments. 

In order to extract the fractional revivals one has to calculate the autocorrelation 
function in a Poincare map for classical cyclotron time differences. Figure displays the 
fractional revivals for several t = {n + p/q) where n is an integer and p/q represents 
an irreducible fraction. Suitable fractions for different fractional revivals are found by 
slicing the circle depicted in figure [7](a). For a full revival the autocorrelation function 
for p/q = needs to be studied, whereas from p/q = 1/2 mirror revivals are identified. 
The next higher fractional revivals with a triangular shape have the time shifts p/q = 0, 
p/q = 1/3, p/q = 2/3 and p/q = 1/2, p/q = 1/6, p/q = 5/6. The absolute value of 
the autocorrelation for these times is shown in figure [3(b) and indicates the times when 
revivals occur. A full revival takes place if the absolute value for p/q = becomes unity 
again, which is the case for t = Tj-^^- The revival for t = l/2T^cv has a minimum for 
p/q = and a maximum for p/q =1/2, explaining the discrepancy between the classical 
centre-of-mass and the quantum mechanical result. Fractional revivals with two peaks 
are visible at times t = l/4Trev and t = S/ATj-ev In this case the absolute value of the 
autocorrelation function has to reach almost y^l/2. Fractional revivals with a triangular 
shape can occur in two different ways. For the first one a maximum is located at the 
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opposite side of the initial starting point. In this case the autocorrelation function 



for p/q = 1/2, p/q = 1/6, p/q = 5/6 is yl/S as shown in the numerical data. As 
second possibility for a triangular revival the same result occurs for p/q = 0, p/q = 1/3, 
p/q = 2/3. Also the next higher revivals can be seen in the numerical data but are not 
presented here to keep clarity. 

3. 5. Effect of impurities 

For graphene based devices, a good understanding of the effect of impurities on the 
electronic motion is required. Suspended graphene has a high mobility leading to 
mean free path lengths of about 2 /im [39]. However, the mobility is reduced by the 
formation of intrinsic ripples which stabilise the two-dimensional system [ID] . The lattice 
deformations involved can be added to the Dirac Hamiltonian by a position dependent 
gauge field [ITl H2] . which can be handled by the propagation algorithm described in 
section [21 Additionally the curvature leads to a spin orbit interaction [13]. Here, we 
study the effect of a scalar potential which is commonly used to model impurities |44j . 
We use a simple model with a Gaussian distributed random potential and analyse 
the propagation of a cyclotron wave-packet through this perturbed environment. The 
potential was generated by the convolution of a random potential t/„ at the grid points 
Vn of a square lattice: 



The potential has a vanishing average expectation value {Uimp{r)) = 0, a variance of 



(t/imp(r)2) = Ul and is correlated by (?7i^p(r)f/i^p(r^)) = ^7^ exp (-|r - 7^17(2^2)) 



with the correlation length ^. A similar potential has already been used in conductance 
calculations for graphene on a honeycomb lattice [S]. We focus on cases where the 
correlation length is longer than the grid spacing leading to a potential which is a 
smooth function with respect to the graphene unit cell. In the following, we study several 
impurity configurations and calculate the propagation of a wave-packet which has the 
same initial shape as a Landau level wave-packet of the clean system. As first observable 
we investigate the autocorrelation function. All the calculated results have in common 
that the autocorrelation function vanishes much faster compared to the clean system. 
In figure [H](a), the grey line shows the autocorrelation function of the clean system with 
prominent revivals whereas the autocorrelation function in the presence of an impurity 
potential decays after a few picoseconds. This phenomenon happens for many different 
setups since the wave-packet simply drifts away. Some impurity configurations lead to 
a wave-packet motion, where the wave-packet returns to its initial position after some 
time. The return results in an increased autocorrelation function but differs from the 
revivals occurring in the clean systems, which refiect the initial occupation of the Landau 





(32) 



levels. 
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Figure 8. (a) Comparison between the autocorrelation function of two different 
impurity configurations (orange and green line, ^ = 42 nm, Uq = 5 meV) and the 
autocorrelation function of the clean system (grey line). The numerical data shows 
that the potential strongly suppresses \C{t) \ of a cyclotron wave-packet (n = 40, a = 1) 
after a few orbits in a magnetic field of -B = 10 T. Revivals in the green line are not 
due to the phase relations between the addends of equation ([27)) but stem from a 
global motion of the wave-packet which returns by chance to the initial position, (b) 
The fidelity M (t) averaged over 10 different impurity configurations for the same initial 
cyclotron wave-packet (n = 40, cr = 1, B = 10 T) feature a similar decay. 



As a measure for the deviation of a disordered system from the clean system we 
calculate the fidelity 

M{t) = |(V;(0)|e'^'*/V^*/^|V^(0))|^ (33) 

We track the overlap of the wave-packet propagated by the unperturbed Hamiltonian H 
and the perturbed Hamiltonian H' = Ti + f/imp- For the free electron gas, expressions are 
known which describe the fidelity decay in disordered systems for differently correlated 
impurity potentials [15]. We obtain similar numerical results, shown in figure [Uj^b). 
For short times the fidelity is approximately unity and then decays exponentially. 
After longer times, the fidelity does not approach zero since there is a non-vanishing 
probability to return back to the initial position. However we cannot infer from 
the autocorrelation function whether wave-packet revivals occur or not, since the 
exponential fidelity decay removes all relevant information from the calculated data. 

Consequently we have to find other observables which signify revivals even if an 
impurity potential is present. One particularly suitable observable is the centre-of-mass 
motion of the wave-packet. In figure [9t^a-c) the centre-of-mass motion of a Dirac wave- 
packet is shown for different impurity potentials. In figure [9](a) the wave-packet does not 
drift away and thus the revivals lie on top of each other. This leads to a revival in the 
autocorrelation function which has nothing in common with the wave-packet revivals 
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Figure 9. (a-c) Centre-of-mass motion of a cyclotron wave- packet {n = 40, a = 1, 
B ~ 10 T) for different impurity configurations = 42 nm, Uq = 5meV). The time- 
evolution shows the collapse of the localised wave-packet by a shrinking spiral. The 
evenly distributed circular wave-packet shows much smother movement leading to a 
connecting line between two revivals, (d) Velocity of the centre-of-mass for the depicted 
wave-packet trajectories from above. The velocity approaches the "speed of light" of 
graphene for each revival. Revivals appear at the same times as in the clean system. 
Although the impurity potential reduces the velocity of the wave-packet for l/2Trev 
and Trev the effect of revivals is still very pronounced in the time-evolution of the 
centre-of-mass velocity. 



of the clean system since the revival time depends crucially on the exact choice of the 
impurity potential. However the majority of the disorder configurations lead to a drift 
of the wave-packet as shown in figure [9]|^b,c). In all cases the wave-packet moves at 
the beginning on a shrinking spiral since the centre-of-mass is located at the position 
with the highest probability density (figure [TOT a)). In this state the centre of the spiral 
already moves. After a few picoseconds the centre-of-mass collapses to the centre of 
a rotating circular wave-packet as shown in figure fTOT b.c). In the following the wave- 
packet is moving without oscillations on a smooth line because the density is equally 
distributed along a circular orbit. As in the clean system, the first revival occurs in the 
form of a mirror revival at time l/2Ti.ev Then the centre-of-mass collapses again and 
produces a full revival at Tj-ev A way to remove the structure of the impurity potential 
is extracting the corresponding velocity of the wave-packet out of the centre-of-mass 
motion 

v=^^{m\r\m), (34) 

which is shown in figure [9](d). The initial wave-packet moves on the circular orbit with 
a velocity slightly smaller than the "speed of light" of graphene. When the centre-of- 
mass collapses this velocity drops down to the drift velocity of the wave-packet which 
is defined by the potential gradient and at least two orders of magnitude smaller. The 
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Figure 10. Revivals and fractional revivals of a cyclotron wave-packet in an impurity 
potential = 42 nm, Uq = 5meV). The red pin marks the quantum mechanical 
centre-of-mass. (a) The initial wave-packet with a Gaussian eigenstate contribution 
n = 40,(T=l. i3 = 10T leads to a Poincare cyclotron time of tO = Tc\. During 
the time evolution the centre-of-mass collapses to the centre of a circular wavefunction 
featuring pronounced revivals. In (b) a quarter revival occurs at t = \/2nTc\ and in (c) 
a half revival aX t = nTd- After a certain time a localised wave-packet emerges again 
and (d) shows a mirror revival for t ~ 2n Td- In this case some part of the wavefunction 
is detached from the localised cyclotron wave-packet, (e) The centre-of-mass collapses 
again and forms fractional revivals, (f) The full revival where the semiclassical and 
the quantum mechanical centre-of-mass coincide occurs at t = AhTd- After this time 
the wavefunction shows pronounced branches. However, the probability distribution 
is clumped around the centre-of-mass showing the characteristics of a revival in the 
velocity expectation value. 

revivals occur at the same time for all the impurity configurations. The substructure of 
the revivals cannot be seen in the centre-of-mass motion of the wave-packet. However 
the full wavefunction shows fractional revivals (figure [TOl). 

In reference [16] the centre-of-mass propagation was used to establish a dipole 
moment. In our case this dipole moment is maximal if the wave-packet revival takes 
place and almost zero if the wave-packet is collapsed similar to the velocity of the 
wave-packet in figure [9]^d). If such a dipole moment is experimentally measurable, the 
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wave-packet revivals can be observed even if disorder displaces the wave-packet far away 
from its initial position. 

4. Conclusion 

In this work we presented a new algorithm for the propagation of wave-packets on a 
sheet of graphene which allows us to calculate the time-evolution in arbitrary shaped 
potentials and inhomogeneous magnetic fields. We studied time-dependent effects like 
the zitterbewegung, which is very pronounced in graphene. We focussed on wave-packet 
revivals, which occur in graphene due to the special structure of the Landau levels. 
We gave a detailed explanation of the centre-of-mass motion and the autocorrelation 
function. The applied semiclassical description quantitatively revealed the collapse of 
the centre of mass as well as the revivals. From the autocorrelation function we extracted 
several different timescales which all signify various revivals. A detailed analysis of the 
autocorrelation function for different Poincare times confirmed fractional revivals which 
have also been seen in our numerical calculations which supported the results throughout 
the whole work. Finally we analysed the effect of impurities on the wave-packet revivals 
by a Gaussian correlated model potential. Our calculations showed that revivals in the 
autocorrelation are strongly suppressed due to the fidelity decay. In contrast to that the 
centre-of-mass revivals are still present if the correlation length of the impurity potential 
is long enough with respect to the average cyclotron radius. The here demonstrated long- 
time accuracy of the propagation algorithm opens the window towards realistic device 
simulations on the micrometre range, also including time-dependent external fields. 
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Appendix A. Classical Propagation 

In the following we study the classical time-evolution of a wave-packet which is located 
around one of the Dirac points and obtain the semiclassical quantisation. In order to be 
localised in position and momentum space, the initial wave-packet ip^r^t) has to have 
a small extent in both coordinates. The expectation value of the position operator and 
the kinetic momentum are defined by 

r{t) = {m\r\m), m = {miMm)- (a.i) 

Niu et. al. have developed a framework to study the time-evolution of those observables 
for arbitrary Hamilton operators [17] which we will apply in the following. As a first 
step we analyse the free solutions of the Hamiltonian and classify the resulting bands 
by the index A. For the Dirac Hamiltonian of graphene we have electron-like (A = +1) 
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and hole-like states (A = —1) with the corresponding energy dispersion Sx = X v h\k\. 
Then it is possible to construct a wave-packet which is restricted to one band. The 
time-evolution of the momentum and the centre of this wave-packet are defined by two 
coupled equations of motion 

, =_i_j=xn.W. (A.2) 

hk= - eE{r) -erx B, (A.3) 

where B = Bbz stands for an isotropic magnetic field perpendicular to the plane and 
E{r) is a position dependent electric potential. The Berry curvature O^x^k) enters the 
equations of motion in momentum space but is zero in the case of a gap-less perfect sheet 
of graphene. As a result of the linear dispersion, relation ( 1A.2I) reduces to r = \vk/\k\ 
and therefore describes a centre-of-mass which always travels with the speed v parallel 
or antiparallel to the wave vector k. For a vanishing electric field the two coupled 
equations of motion can be solved. The momentum of the wave-packet is controlled by 
the differential equation 

k 

hk= -ev— X B(r), (A.4) 
\k\ 

given by (1A.3I) which is solved by 

y sm[uJct + ifi) J 

where the initial momentum of the wave-packet is given by ko in direction if. The 
cyclotron frequency depends on the momentum and is given by 

uj, = —-. A.6 

hko 

The resulting real space propagation flA.2p is given by 

, , / sinftUct 0?) \ / A 

r=ro + A/J \' . A.7 

Y — cos(co'ct: + ip) J 

Thus electron-like and hole-like states propagate in opposite directions on orbits with 
the radius 

4 ^ (A.8, 

In order to obtain the quantisation of the Landau levels from the closed cyclotron orbits 
we have to study the phase changes along the trajectories. The main contribution 
is given by the classical action known from the electron case, while an additional 
contribution comes from Berry's phase. Although the Berry curvature is vanishing. 
Berry's phase [IH] is still present and defined by the vector potential 

Ax{k) = ~i{uxik)\'^kUx{k)), (A.9) 

and the spinor parts ux{k) of the free solutions of the Hamiltonian. The free solutions 
are given by 

u^{k) = (I) , (A.IO) 
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where (p is the direction of the vector k. Hence we get a vector potential of 

which is used to calculate the additional phase-change 

Tx{uj) = <f Ax{k)dk (A.12) 

J UJ 

for an arbitrary adiabatic transition on the path uj{t) in momentum space. The result 
is Berry's phase 

r±(cu) = ±71 (A.13) 

for a single closed orbit on graphene [2]. The adiabatic transition in momentum space 
is automatically fulfilled in a magnetic field, since all trajectories are cyclotron orbits 
and thus sufficiently smooth. The classical phase change and Berry's phase have to be 
added up to build a new EBK quantisation condition which reads [19] 



Ufkx dk=2n'4- fn' + ^ - (A.14) 

where v is the Maslov index depending on the caustics traversed on one cyclotron 
orbit (here u = 2). Plugging the formula for the cyclotron motion in the quantisation 
condition the resulting Landau levels are given by 

nkl = 27x—[n' + (A.15) 

If we define the Landau level index n differently for hole-like and electron-like states the 
momentum quantisation yields 



/ eB 

kf) = sign{n)\\2—\n\. (A. 16) 

All the other quantised parameters follow from the relations calculated before 



Sn = V sign(n)A/2ei?/i|n|, (A. 17) 



sign(n)W2— |n|, (A.18) 



/ eB 

Un = v sign(n) J — — - if n^O. (A. 19) 

V 2n\n\ 
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